classdef getDensity < handle
	%GETDENSITY computes bootstrapped density estimates and full stats for two
	%groups. You can pass a set of cases to the function and it will compute
	%statistics automagically for all cases as well. There are numerous
	%settings:
	% x - our first data group, plotted along the abscissa
	% y - our second data group, plotted along the ordinate
	% nboot - number of bootstrap iterations default: 1000
	% avgFunction - function handle to bootstrap default: @mean
	% alpha - P value to use for all statistics
	% legendtxt - A cell list of the main X and Y data, i.e. {'Control','Drug'}
	% dogauss - Add gaussian fits to the histograms? default: true
	% columnlabels - if X/Y have more than 1 column, you can detail column names here
	% dooutlier - remove outliers using 'grubb', 'rosner', 3SD' or 'none'
	% addjitter - add jitter to scatterplot? 'x', 'y', 'both', 'equal', 'none'
	% cases - a celllist of case names the same length as X/Y
	% jitterfactor - how much to scale the jitter, larger values = less jitter
	% showoriginalscatter - show unjittered data too?
	
	properties
		%> comments to add to this data comparison
		comments = {''}
		%> our first data group, plotted along the abscissa
		x = []
		%> our second data group, plotted along the ordinate
		y = []
        %> a cell of names for each row
		xrownames = []
		%> a cell of names for each row
		yrownames = []
		%> number of bootstrap iterations default: 1000
		nboot = 1000
		%> function handle to bootstrap default: @mean
		avgFunction = @mean
		%> P value to use for all statistics
		alpha = 0.05
		%> A cell list of the main X and Y data, i.e. {'Control','Drug'}
		legendtxt = {'Group1','Group2'}
		%> if X/Y have more than 1 column, you can detail column names here
		columnlabels = {'DataSet1'}
        %> Add gaussian fits to the histograms? default: true
		dogauss = true
		%> add jitter to scatterplot? 'x', 'y', 'both', 'equal', 'none'
		addjitter = 'none'
		%>  celllist of case names the same length as X/Y
		cases = []
		%> are our cases nominal or ordinal?
		nominalcases = true
		%> remove outliers using 'grubb', 'rosner', 'limit', '3SD' or 'none'
		dooutlier = 'none'
		%> rosner outlier number of outliers
		rosnern = 2
		%> p to use for outlier removal for rosner and grubb
		outlierp = 0.001
		%> outlier limit for limit method
		outlierlimit = Inf
		%> scalefactor for jitter of scatter plot
		jitterfactor = 100
		%> do we spit out each plot in its own figure for easy export?
		singleplots = false
		%> do we show the original points in the scatter before scattering them?
		showoriginalscatter = false
		%> if we are passed data on construction, run straight away
		autorun = false
		%> how to plot the histogram
		bartype = 'grouped'
		%> width of histogram bars
		barwidth = 1.75
		%> are we verbose printing to command line?
		verbose = true
		%> do we scale the axes to be equal for the scatterplots
		scaleaxes = true
		%> do we show the unity line for the scatterplot?
		showunityline = true
		%>number of bins for histogram
		nBins = 10
		%> centre histogram bins
		centreBins = true
		%>density kernel to use, 'normal', 'box', 'triangle', 'epanechnikov'
		densityKernel = 'normal'
		%> density function, 'pdf', 'cdf','icdf', 'survivor', 'cumhazard'
		densityFunction = 'pdf'
		%> density bounds; 'unbounded', 'positive' or 2-element vector for lower and upper
		densityBounds = 'unbounded'
		%> for multiple columns only run a subset; empty runs all.
		index = []
		%> normalise (Z-score) values to plot on scatter?
		normaliseScatter = false
		%> show distribution plot underneath box plot?
		doDistribution = true
        %> do a median CI, takes longer if dataset is large...
        doMedianCI = true
	end
	
	properties (Dependent = true, SetAccess = private, GetAccess = public)
		%> is the x and y data groups equal length, suggestive of paired data?
		isDataEqualLength = true;
		%> number of columns
		nColumns = []
		%> the unique cases
		uniquecases = []
	end
	
	properties (Dependent = true, SetAccess = private, GetAccess = private)
		%> the layout when plotting mutliple plots
		plotX = 3
		%> layout when plotting mutliple plots
		plotY = 2
	end
	
	properties (SetAccess = private, GetAccess = public)
		%> the structure returned from run, saved here
		runStructure
		%> autogenerated dataset for first group
		xdata = []
		%> autogenerated dataset for second group
		ydata = []
	end
	
	properties (SetAccess = private, GetAccess = private)
		%> our panel object, we use panel instead of subplot as it is more powerful
		pn
		workindex
		allowedProperties = []
		ignoreProperties = '(isDataEqualLength|nColumns|uniquecases|runStructure|h|pn|plotX|plotY)';
	end
	
	events (Hidden = false, ListenAccess = private, NotifyAccess = private)
		checkData %trigger this event to perform checks on data and cases etc.
	end
	
	%=======================================================================
	methods %------------------PUBLIC METHODS
		%=======================================================================
		
		% ===================================================================
		%> @brief Class constructor
		%>
		%> More detailed description of what the constructor does.
		%>
		%> @param args are passed as a structure of properties which is
		%> parsed.
		%> @return instance of class.
		% ===================================================================
		function obj = getDensity(varargin)
			
			addlistener(obj,'checkData',@obj.doCheckData); %use an event to keep data accurate
			
			if nargin>0
				obj.parseArgs(varargin, obj.allowedProperties);
			end
			
			if obj.autorun == true && obj.haveData == true
				obj.run;
			end
		end
		
		% ===================================================================
		%> @brief Our main RUN method
		%>
		%> This performs the analysis and plots the results for the stored data.
		%>
		%> @param obj the class object automatically passed in
		%> @return outs a structure with the analysis results if needed.
		% ===================================================================
		function outs = run(obj)
			notify(obj,'checkData');
			if isempty(obj.x)
				error('You haven''t supplied any data yet!')
			end
			warning('off', 'stats:lillietest:OutOfRangePLow')
			warning('off', 'stats:lillietest:OutOfRangePHigh')
			warning('off', 'stats:jbtest:PTooSmall')
			warning('off', 'stats:jbtest:PTooBig')
			if ~isempty(obj.index)
				nVals = length(obj.index);
				obj.workindex = obj.index;
			else
				nVals = size(obj.x,2);
				obj.workindex = 1:nVals;
			end
			for i = 1:nVals
				idx=obj.workindex(i);
				t0 = tic;
				xcol=obj.x(:,idx);
				ycol=obj.y(:,idx);
				c1 = [0.2 0.2 0.2];
				c2 = [0.8 0.2 0.2];
				casesLocal = obj.cases;
				
				outs.x = obj.x;
				outs.y = obj.y;
				outs.cases = casesLocal;
				
				xouttext='';
				youttext='';
				
				fieldn = obj.columnlabels{idx};
				fieldn = regexprep(fieldn,'\s+','_');
				
				h=figure;
				outs.(fieldn).h = h;
				set(h,'Color',[1 1 1]);
				t = [obj.columnlabels{:}];
				set(h,'name',t);
				
				pn = panel(h);
				if obj.isDataEqualLength == false
					if exist('figpos','file'); figpos(1,[1000 1000]); end
					delete(gca)
					pn.pack('v',[0.5 0.5]);
					pn(1).pack('h',[1/3 1/3 -1])
					pn(2).pack('h',[1/3 -1])
				else
					if exist('figpos','file'); figpos(1,[1200 1000]); end
					delete(gca)
					pn.pack('v',[0.5 0.5]);
					pn(1).pack('h',[1/3 1/3 -1])
					pn(2).pack('h',[1/3 1/3 -1])
				end
				
				if any(isnan(xcol)) %remove any nans
					xcol(isnan(xcol))=[];
				end
				if any(isnan(ycol))
					ycol(isnan(ycol))=[];
				end
				
				[xcol,ycol,casesLocal,xouttext,youttext] = obj.removeOutliers(xcol,ycol,casesLocal,xouttext,youttext);
				
				xmean=mean(xcol, 'omitnan');
				xmedian=median(xcol, 'omitnan');
				ymean=mean(ycol, 'omitnan');
				ymedian=median(ycol, 'omitnan');
				xstd=std(xcol, 'omitnan');
				ystd=std(ycol, 'omitnan');
				xstderr=obj.stderr(xcol,'SE',1);
				ystderr=obj.stderr(ycol,'SE',1);
				
				%==========================================DO HISTOGRAM
				px = 1;
				py = 1;
				pn(py,px).select();
				
				if ystd > 0
					[n1,b1]=histcounts(xcol);
					[n2,b2]=histcounts(ycol);
					bw = min([diff(b2) diff(b1)]);
					bl = [min([min(xcol) min(ycol)]) max([max(xcol) max(ycol)])];
					outs.(fieldn).histo1 = histogram(xcol,10,'BinWidth',bw,...
						'BinLimits',bl,'EdgeColor','none','FaceAlpha',0.6);
					outs.(fieldn).histo1.FaceColor = c1;
					xn = outs.(fieldn).histo1.Values;
					hbinsx = outs.(fieldn).histo1.BinEdges;
					hold on
					outs.(fieldn).histo2 = histogram(ycol,10,'BinWidth',bw,...
						'BinLimits',bl,'EdgeColor','none','FaceAlpha',0.6);
					outs.(fieldn).histo2.FaceColor = c2;
					yn = outs.(fieldn).histo2.Values;
					hbinsy = outs.(fieldn).histo2.BinEdges;
				else
					outs.(fieldn).histo1 = histogram(xcol,10,'EdgeColor','none','FaceAlpha',0.6);
					outs.(fieldn).histo1.FaceColor = c1;
					xn = outs.(fieldn).histo1.Values;
					hbinsx = outs.(fieldn).histo1.BinEdges;
				end
				
				lim=ylim;
				text(xmean,lim(2),'\downarrow','Color',c1,'HorizontalAlignment','center',...
					'VerticalAlignment','bottom','FontSize',15,'FontWeight','bold');
				text(xmedian,lim(2),'\nabla','Color',c1,'HorizontalAlignment','center',...
					'VerticalAlignment','bottom','FontSize',15,'FontWeight','bold');
				if ystd > 0
					text(ymean,lim(2),'\downarrow','Color',c2,'HorizontalAlignment','center',...
						'VerticalAlignment','bottom','FontSize',15,'FontWeight','bold');
					text(ymedian,lim(2),'\nabla','Color',c2,'HorizontalAlignment','center',...
						'VerticalAlignment','bottom','FontSize',15,'FontWeight','bold');
				end
				if obj.dogauss == true
					hold on
					if length(find(xn>0))>1
						try
							f1=fit(hbinsx(1:end-1)',xn','gauss1');
							plot(f1,'Color',c1);
						catch
							disp('Couldn''t fit gaussian to histogram 1')
						end
					end
					if ystd > 0 && length(find(yn>0))>1
						try
							f2=fit(hbinsy(1:end-1)',yn','gauss1');
							plot(f2,'Color',c2);
						catch
							disp('Couldn''t fit gaussian to histogram 2')
						end
					end
					legend off
					hold off
				end
				
				grid on; box on;
				pn(py,px).xlabel(obj.columnlabels{idx});
				pn(py,px).ylabel('Number of Units');
				pn(py,px).title('Histogram');
                drawnow;
				
				%==========================================DO BOX PLOTS
				px = 2;
				py = 1;
				pn(py,px).select()
				hold on
				if exist('distributionPlot','file') && obj.doDistribution
					distributionPlot({xcol,ycol},'colormap',parula,'showMM',2,...
						'xValues',[0.625 1.625],'histOri','center','distWidth',0.25);
				elseif obj.doDistribution
					obj.violin(0.5,xcol,'withmdn',false,'facecolor',c1,'side','both','scaling',0.4);
					obj.violin(1.5,ycol,'withmdn',false,'facecolor',c2,'side','both','scaling',0.4);
				end
				if obj.isDataEqualLength
					boxplot([xcol ycol],'positions',[1 2],'notch','on','whisker',1.5,...
						'labels',obj.legendtxt,'colors','k',...
						'boxstyle','outline','medianstyle','line',...
						'widths',0.5,'symbol','k.','OutlierSize',3,'Jitter',0.25);
				else
					hold on
					boxplot(xcol,'positions',1,'notch','on','whisker',1.5,...
						'labels',obj.legendtxt{1},'colors',c1,...
						'boxstyle','outline','medianstyle','line',...
						'widths',0.5,'symbol','k.','OutlierSize',2,'Jitter',0.25);
					boxplot(ycol,'positions',2,'notch','on','whisker',1.5,...
						'labels',obj.legendtxt{2},'colors',c2,...
						'boxstyle','outline','medianstyle','line',...
						'widths',0.5,'symbol','k.','OutlierSize',2,'Jitter',0.25);
				end
				pn(py,px).ylabel(obj.columnlabels{idx});
				axis tight
				if obj.doDistribution
					xlim([0.1 2.4]);
				else
					xlim([0.5 2.5]);
				end
				set(gca,'XTick', [1 2],'XTickLabel', obj.legendtxt)
				pn(py,px).title('Box / Density Plots')
				hold off
				grid on; box on
				yl1 = ylim; %we check this against the next plot below
            drawnow;
				
				%==========================================DO SCATBOX PLOTS
				px = 3;
				py = 1;
				pn(py,px).select()
				
				if obj.isDataEqualLength
					% 					if ~isempty(obj.uniquecases) && ~isempty(casesLocal)
					% 						for jj = 1:length(obj.uniquecases)
					% 							caseidx = ismember(casesLocal,obj.uniquecases{jj});
					% 							xtmp(:,jj) = xcol(caseidx);
					% 							ytmp(:,jj) = ycol(caseidx);
					% 						end
					% 						notBoxPlot([xtmp ytmp],[1 1 2 2]);
					% 					else
					% 						notBoxPlot([xcol ycol],[1 1]);
					% 					end
					obj.notBoxPlot([xcol ycol],[1 2]);
					set(gca,'XTick', [1 2],'XTickLabel', obj.legendtxt)
					pn(py,px).ylabel(obj.columnlabels{idx});
					
				else
					
					obj.notBoxPlot(xcol,1);
					hold on
					obj.notBoxPlot(ycol,2);
					set(gca,'XTick', [1 2],'XTickLabel', obj.legendtxt)
					pn(py,px).ylabel(obj.columnlabels{idx});
					
				end
				xlim([0 3]);
				pn(py,px).title('ScatterBox Plots')
				hold off
				grid on; box on
				yl2 = ylim;
				
				%==========================================EQUALISE Y AXIS
				ym1 = min([yl1(1) yl2(1)]);
				ym2 = max([yl1(2) yl2(2)]);
				pn(1,2).select()
				ylim([ym1 ym2]);
				pn(1,3).select()
				ylim([ym1 ym2]);
                drawnow;
				
				%==========================================DO Correlation SCATTER
				xcolout = xcol;
				ycolout = ycol;
				if ystd > 0 && obj.isDataEqualLength
					px = 1;
					py = 2;
					pn(py,px).select()
					[r,p]=corr(xcol,ycol);
					[r2,p2]=corr(xcol,ycol,'type','spearman');
					
					xrange = max(xcol) - min(xcol);
					yrange = max(ycol) - min(ycol);
					range = max(xrange,yrange);
					
					
					if obj.addjitter == true
						obj.addjitter = 'both';
					end
					switch obj.addjitter
						case 'x'
							sc = true;
							xcolout = obj.jitterData(xcol);
							ycolout = ycol;
						case 'y'
							sc = true;
							xcolout = xcol;
							ycolout = obj.jitterData(ycol);
						case 'equal'
							sc = true;
							[xcolout,jitter] = obj.jitterData([xcol ycol]);
							ycolout = obj.jitterData([ycol xcol],jitter);
						case 'both'
							sc = true;
							xcolout = obj.jitterData(xcol);
							ycolout = obj.jitterData(ycol);
						otherwise
							sc = false;
							xcolout = xcol;
							ycolout = ycol;
					end
					if obj.normaliseScatter == true
						xcolout = zscore(xcolout);
						ycolout = zscore(ycolout);
						mn = min(min(zscore(xcol)),min(zscore(ycol)));
						mx = max(max(zscore(xcol)),max(zscore(ycol)));
					else
						mn = min(min(xcol),min(ycol));
						mx = max(max(xcol),max(ycol));
					end
					if obj.scaleaxes == true
						axrange = [(mn - (mn/20)) (mx + (mx/20)) (mn - (mn/20)) (mx + (mx/20))];
					else
						axrange = [min(xcol) min(xcol)+range min(ycol) min(ycol)+range];
					end
					if ~isempty(casesLocal)
						t = 'Group: ';
						for jj = 1:length(obj.uniquecases)
							t = [t num2str(jj) '=' obj.uniquecases{jj} ' '];
						end
					else
						colours = [0 0 0];
						t = '';
					end
					
					hold on
					if ~isempty(casesLocal)
						gscatter(xcolout,ycolout,casesLocal,'krbgmyc','o');
					else
						scatter(xcolout,ycolout,repmat(80,length(xcolout),1),'o','MarkerEdgeColor', [0 0 0],...
							'MarkerFaceColor', [0.65 0.65 0.65],'MarkerFaceAlpha',0.2,'MarkerEdgeAlpha',0.75);
					end
					try %#ok<TRYNC>
						h = lsline;
						set(h,'LineStyle','--','LineWidth',2)
						if obj.showunityline == true
							if abs(r2) < 0.1
								h = line([mn mx],[mn mx]);
								set(h,'Color',[0.7 0.7 0.7],'LineStyle','-.','LineWidth',2)
								h = line([mx mn],[mn mx]);
								set(h,'Color',[0.7 0.7 0.7],'LineStyle','-.','LineWidth',2)
							elseif r2 >= 0
								h = line([mn mx],[mn mx]);
								set(h,'Color',[0.7 0.7 0.7],'LineStyle','-.','LineWidth',2)
							else
								h = line([mx mn],[mn mx]);
								set(h,'Color',[0.7 0.7 0.7],'LineStyle','-.','LineWidth',2)
							end
						end
					end
					if sc == true && obj.showoriginalscatter == true
						scatter(xcol,ycol,repmat(80,length(xcol),1),'k.','MarkerFaceColor',[0.5 0.5 0.5],'MarkerEdgeColor',[0.5 0.5 0.5]);
					end
					axis square
					axis(axrange);
					if obj.scaleaxes == true
						set(gca,'XTick',get(gca,'YTick'),'XTickLabel',get(gca,'YTickLabel'));
					end
					if obj.normaliseScatter == true
						pn(py,px).xlabel(['ZScore ' obj.legendtxt{1}])
						pn(py,px).ylabel(['ZScore ' obj.legendtxt{2}])
					else
						pn(py,px).xlabel(obj.legendtxt{1})
						pn(py,px).ylabel(obj.legendtxt{2})
					end
					pn(py,px).title(['Prson:' sprintf('%0.2g',r) '(p=' sprintf('%0.4g',p) ') | Spman:' sprintf('%0.2g',r2) '(p=' sprintf('%0.4g',p2) ') ' t]);
					hold off
					grid on; box on
					set(gca,'Layer','top');
                    drawnow;
				end
				
				%============================Lets measure statistics
				t=['Mn/Mdn: ' sprintf('%0.3g', xmean) '\pm' sprintf('%0.3g', xstderr) '/' sprintf('%0.3g', xmedian) ' | ' sprintf('%0.3g', ymean) '\pm' sprintf('%0.3g', ystderr) ' / ' sprintf('%0.3g', ymedian)];
				
				if ystd > 0
					if length(xcol) == length(ycol)
						[h,p1]=ttest2(xcol,ycol,obj.alpha);
					else
						[h,p1]=ttest2(xcol,ycol,obj.alpha);
					end
					[p2,h]=ranksum(xcol,ycol,'alpha',obj.alpha);
					if obj.isDataEqualLength
						[h,p3]=ttest(xcol,ycol,obj.alpha);
						[p4,h]=signrank(xcol,ycol,'alpha',obj.alpha);
						[p5,h]=signtest(xcol,ycol,'alpha',obj.alpha);
					else
						p3 = [];
						p4 = [];
						p5 = [];
					end
					[h,p6]=jbtest(xcol,obj.alpha);
					[h,p7]=jbtest(ycol,obj.alpha);
					if length(xcol)>4
						[h,p8]=lillietest(xcol);
					else
						p8=NaN;
					end
					if length(ycol)>4
						[h,p9]=lillietest(ycol);
					else
						p9=NaN;
					end
					[h,p10]=kstest2(xcol,ycol,obj.alpha);
				else
					[h,p1]=ttest(xcol,obj.alpha);
					[p2,h]=signrank(xcol,0,'alpha',obj.alpha);
					p3=[];
					p4=[];
					[p5,h]=signtest(xcol,0,'alpha',obj.alpha);
					[h,p6]=jbtest(xcol,obj.alpha);
					p7=0;
					if length(xcol)>4
						[h,p8]=lillietest(xcol);
					else
						p8=NaN;
					end
					p9=[];
					p10=kstest(xcol);
				end
				
				outs.(fieldn).ttest = p1;
				outs.(fieldn).ranksum = p2;
				outs.(fieldn).ttest2 = p3;
				outs.(fieldn).signrank = p4;
				outs.(fieldn).signtest = p5;
				outs.(fieldn).jbtestx = p6;
				outs.(fieldn).jbtesty = p7;
				outs.(fieldn).lillietestx = p8;
				outs.(fieldn).lillietesty = p9;
				outs.(fieldn).kstest = p10;
				
				t=[t '\newlineT-test: ' sprintf('%0.3g', p1) '\newline'];
				t=[t 'Wilcox: ' sprintf('%0.3g', p2) '\newline'];
				if exist('p3','var') && ystd > 0
					t=[t 'Pair ttest: ' sprintf('%0.3g', p3) '\newline'];
					t=[t 'Pair wilcox: ' sprintf('%0.3g', p4) '\newline'];
					t=[t 'Pair sign: ' sprintf('%0.3g', p5) '\newline'];
				elseif exist('p5','var')
					t=[t 'Sign: ' sprintf('%0.3g', p5) '\newline'];
				end
				t=[t 'Jarque-Bera: ' sprintf('%0.3g', p6) ' / ' sprintf('%0.3g', p7) '\newline'];
				t=[t 'Lilliefors: ' sprintf('%0.3g', p8) ' / ' sprintf('%0.3g', p9) '\newline'];
				t=[t 'KSTest: ' sprintf('%0.3g', p10) '\newline'];
				
                disp('Calculating bootstrap CI...')
                options = statset('UseParallel',true);
                t1 = tic;
				[xci,xpop]=bootci(obj.nboot,{obj.avgFunction,xcol},'alpha',obj.alpha,'Options',options);
				[yci,ypop]=bootci(obj.nboot,{obj.avgFunction,ycol},'alpha',obj.alpha,'Options',options);
                fprintf('---> took: %.2g seconds\n',toc(t1));
                disp('Calculating bootstrap mean...')
                t1b = tic;
				xmean = obj.avgFunction(xpop);
				ymean = obj.avgFunction(ypop);
                fprintf('---> took: %.2g seconds\n',toc(t1b));
                
                t=[t 'BootStrap: ' sprintf('%0.3g', xci(1)) ' < ' sprintf('%0.3g', xmean) ' > ' sprintf('%0.3g', xci(2)) ' | ' sprintf('%0.3g', yci(1)) ' < ' sprintf('%0.3g', ymean) ' > ' sprintf('%0.3g', yci(2))];
				
                if obj.doMedianCI
                    try %#ok<*TRYNC>
                        disp('Calculating median bootstrap CI...')
                        t2 = tic;
                        [xxci,xxpop]=bootci(obj.nboot,{@median,xcol},'alpha',obj.alpha,'Options',options);
                        [yyci,yypop]=bootci(obj.nboot,{@median,ycol},'alpha',obj.alpha,'Options',options);
                        disp('Calculating bootstrap median...')
                        xxmean = median(xxpop);
                        yymean = median(yypop);
                        fprintf('---> took: %.2g seconds\n',toc(t2));
                        t=[t '\newlineMedian BS: ' sprintf('%0.3g', xxci(1)) ' < ' sprintf('%0.3g', xxmean) ' > ' sprintf('%0.3g', xxci(2)) ' | ' sprintf('%0.3g', yyci(1)) ' < ' sprintf('%0.3g', yymean) ' > ' sprintf('%0.3g', yyci(2))];
                    end   
                    try
                        disp('Calculating geomean bootstrap CI...')
                        t3 = tic;
                        [xxxci,xxxpop]=bootci(obj.nboot,{@geomean,xcol},'alpha',obj.alpha,'Options',options);
                        [yyyci,yyypop]=bootci(obj.nboot,{@geomean,ycol},'alpha',obj.alpha,'Options',options);
                        disp('Calculating bootstrap geomean...')
                        xxxmean = geomean(xxxpop);
                        yyymean = geomean(yyypop);
                        fprintf('---> took: %.2g seconds\n',toc(t3));
                        t=[t '\newlineGeomean BS: ' sprintf('%0.3g', xxxci(1)) ' < ' sprintf('%0.3g', xxxmean) ' > ' sprintf('%0.3g', xxxci(2)) ' | ' sprintf('%0.3g', yyyci(1)) ' < ' sprintf('%0.3g', yyymean) ' > ' sprintf('%0.3g', yyyci(2))];
                    end
                end
				
                disp('Calculating KS Density function...')
                t4 = tic;
				[fx,xax]=ksdensity(xpop, 'kernel', obj.densityKernel,...
					'function', obj.densityFunction, 'support', obj.densityBounds);
				[fy,yax]=ksdensity(ypop, 'kernel', obj.densityKernel,...
					'function', obj.densityFunction, 'support', obj.densityBounds);
                fprintf('---> took: %.2g seconds\n',toc(t4));
                
                clear xpop ypop xxpop yypop xxxpop yyypop
				
				%==========================================DO CDF
				if obj.isDataEqualLength
					px = 2;
					py = 2;
				else
					px = 1;
					py = 2;
				end
				pn(py,px).select();
				
				hold on
				[f,x,flo,fup] = ecdf(xcol,'alpha',obj.alpha);
				flo(1)=0;flo(end)=flo(end-1); fup(1)=0; fup(end)=fup(end-1);
				patch('XData',[x;flipud(x)],'YData',[flo;flipud(fup)],'FaceColor',[0 0 0],'FaceAlpha',0.2,'EdgeColor','none');
				stairs(x,f,'k','LineWidth',2);
				[f,x,flo,fup] = ecdf(ycol,'alpha',obj.alpha);
				flo(1)=0;flo(end)=flo(end-1); fup(1)=0; fup(end)=fup(end-1);
				patch('XData',[x;flipud(x)],'YData',[flo;flipud(fup)],'FaceColor',[1 0 0],'FaceAlpha',0.2,'EdgeColor','none');
				stairs(x,f,'r','LineWidth',2);
				hold off
				
				grid on; box on
				pn(py,px).title(['Cumulative Distribution Function, p=' num2str(obj.alpha)]);
				pn(py,px).xlabel(obj.columnlabels{idx});
                drawnow;
				
				%==========================================Do DENSITY
				if obj.isDataEqualLength
					px = 3;
					py = 2;
				else
					px = 2;
					py = 2;
				end
				pn(py,px).select();
				plot(xax,fx,'k-','linewidth',1.5);
				if ystd > 0
					hold on
					plot(yax,fy,'r-','linewidth',1.5);
					hold off
				end
				axis tight
				title(['BootStrap Density Plots; using: ' func2str(obj.avgFunction)]);
				box on
				grid on
				set(gca,'Layer','top');
				
				supt=[obj.columnlabels{idx} ' # = ' num2str(length(xcol)) '[' num2str(length(obj.x(:,i))) '] & ' num2str(length(ycol)) '[' num2str(length(obj.y(:,i))) ']'];
				
				if ~isempty(xouttext) && ~strcmp(xouttext,' ') && length(xouttext)<25
					supt=[supt ' | ' obj.dooutlier ': 1 = ' xouttext];
				end
				if ~isempty(youttext) && ~strcmp(youttext,' ') && length(xouttext)<25
					supt=[supt ' | 2 = ' youttext];
				end
				
				supt = [supt '\newline' obj.comments];
				
				if exist('suptitle','file')
					h = suptitle(supt);
					set(h,'FontName','Menlo','FontSize',12,'FontWeight','bold')
				else
					disp(supt);
				end
				
				xl=xlim;
				yl=ylim;
				
				xfrag=(xl(2)-xl(1))/40;
				yfrag=(yl(2)-yl(1))/40;
				
				h=line([xci(1),xmean,xci(2);xci(1),xmean,xci(2);],[yl(1),yl(1),yl(1);yl(2),yl(2),yl(2)]);
				set(h,'Color',[0.5 0.5 0.5],'LineStyle','--');
				
				h=line([yci(1),ymean,yci(2);yci(1),ymean,yci(2);],[yl(1),yl(1),yl(1);yl(2),yl(2),yl(2)]);
				set(h,'Color',[1 0.5 0.5],'LineStyle',':');
				
				text(xl(1)+xfrag,yl(2)-yfrag,t,'FontSize',10,'FontName','Menlo','FontWeight','bold','VerticalAlignment','top');
				%legend(obj.legendtxt);
				
				pn.select('all')
				pn.fontname = 'Menlo';
				pn.fontsize = 10;
				pn.margin = 20;
				
				outs.(fieldn).pn = pn;
				outs.(fieldn).xcol = xcol;
				outs.(fieldn).ycol = ycol;
				outs.(fieldn).xmean = xmean;
				outs.(fieldn).xmedian = xmedian;
				outs.(fieldn).xstd = xstd;
				outs.(fieldn).xse = xstderr;
				outs.(fieldn).xci = xci;
				outs.(fieldn).ymean = ymean;
				outs.(fieldn).ymedian = ymedian;
				outs.(fieldn).ystd = ystd;
				outs.(fieldn).yse = ystderr;
				outs.(fieldn).yci = yci;
				outs.(fieldn).xcolout = xcolout;
				outs.(fieldn).ycolout = ycolout;
				outs.(fieldn).text = t;
				
				fprintf('\n---> getDensity Total Computation time took: %.2g seconds\n',toc(t0));
				
				obj.pn = pn;
				
				if obj.singleplots == true
					obj.doSinglePlots(pn);
				end
				
				if obj.isDataEqualLength && ~isempty(obj.uniquecases) && ~isempty(casesLocal)
					for jj = 1:length(obj.uniquecases)
						caseidx = ismember(casesLocal,obj.uniquecases{jj});
						xtmp = xcol(caseidx);
						ytmp = ycol(caseidx);
						name = ['Case_' obj.uniquecases{jj}];
						otmp = obj.toStructure();
						otmp.x = xtmp;
						otmp.y = ytmp;
						otmp.cases = cell(0,0);
						otmp.columnlabels{1} = [otmp.columnlabels{i} ' ' name];
						otmp.autorun = false;
						otmp.verbose = false;
						gdtmp = getDensity(otmp);
						outtmp = gdtmp.run;
						outs.(fieldn).(name)=outtmp;
					end
					figure(outs.(fieldn).h);
				end
			end
			obj.runStructure = outs;
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function scatterColumns(obj,c1,c2)
			
			obs1 = obj.xdata.Properties.VarNames{c1};
			obs2 = obj.xdata.Properties.VarNames{c2};
			
			d1 = obj.xdata.(obs1);
			d2 = obj.xdata.(obs2);
			d3 = obj.ydata.(obs1);
			d4 = obj.ydata.(obs2);
			
			if obj.addjitter == true
				obj.addjitter = 'both';
			end
			switch obj.addjitter
				case 'x'
					sc = true;
					d1 = obj.jitterData(d1);
					d3 = obj.jitterData(d3);
				case 'y'
					sc = true;
					d2 = obj.jitterData(d2);
					d4 = obj.jitterData(d4);
				case 'equal'
					sc = true;
					[d1, jitter] = obj.jitterData([d1 d3]);
					d3 = obj.jitterData([d3 d1],jitter);
					[d2, jitter] = obj.jitterData([d2 d4]);
					d4 = obj.jitterData([d4 d2],jitter);
				case 'both'
					sc = true;
					d1 = obj.jitterData(d1);
					d3 = obj.jitterData(d3);
					d2 = obj.jitterData(d2);
					d4 = obj.jitterData(d4);
				otherwise
					sc = false;
			end
			
			h=figure;
			figpos(1,[1000,1000]);
			set(h,'Color',[1 1 1]);
			hold on
			scatter(d1, d2, repmat(80,length(d1),1),[0 0 0],'MarkerFaceColor','none',...
				'DisplayName',obj.legendtxt{1});
			scatter(d3, d4, repmat(80,length(d3),1),[0.8 0 0],'MarkerFaceColor','none',...
				'DisplayName',obj.legendtxt{2});
			set(gca, 'FontSize', 14);
			hold off
			box on
			grid on
			xlabel([obs1])
			ylabel([obs2])
			title(['Scatterbox plot comparing ' obs1 ' against ' obs2 ' -- jitter:' obj.addjitter]);
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function [data,jitter] = jitterData(obj,datain,jitter)
			
			if ~exist('jitter','var')
				jitter = [];
			end
			
			if size(datain,2) > size(datain,1)
				datain = datain';
			end
			
			if size(datain,2) > 1
				
				other = datain(:,2);
				datain = datain(:,1);
				range1 = max(datain) - min(datain);
				range2 = max(other) - min(other);
				if isempty(jitter)
					jitter = (randn(length(datain),1))*(max(range1,range2)/obj.jitterfactor);
				end
			else
				range = max(datain) - min(datain);
				if isempty(jitter)
					jitter = (randn(length(datain),1))*(range/obj.jitterfactor);
				end
			end
			
			data = datain + jitter;
			
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function value = get.nColumns(obj)
			value = size(obj.x,2);
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function value = get.isDataEqualLength(obj)
			if size(obj.x,1) == size(obj.y,1)
				value = true;
			else
				value = false;
			end
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function value = get.plotX(obj)
			if obj.isDataEqualLength
				value = 3;
			else
				value = 3;
			end
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function value = get.plotY(obj)
			if obj.isDataEqualLength
				value = 2;
			else
				value = 2;
			end
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function value = get.uniquecases(obj)
			if ~isempty(obj.cases)
				value = getlabels(obj.cases);
			else
				value = [];
			end
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function set.x(obj,value)
			if isstruct(value) && length(value) == 1
				f = fieldnames(value);
				firstLength = [];
				x = [];
				names = {};
				for i = 1:length(f)
					col = value.(f{i});
					if isnumeric(col) && isvector(col)
                        if any(size(col)==1) && ~iscolumn(col)
							col = col';
                        end
                        
						if isempty(firstLength)
							firstLength = length(col);
						end
						if length(col) == firstLength
							x(:,end+1) = col;
							names{end+1} = f{i};
						end
					end
				end
				if ~isempty(x)
					obj.x = x;
					obj.columnlabels = names;
				end
			elseif isnumeric(value)
                if any(size(value)==1) && ~iscolumn(value)
                    value=value';
                end
				value(isnan(value)) = []; %purge nans
				obj.x = value;
			elseif istable(value)
				obj.x = table2array(value);
			else
				warning('x input data isn''t valid, must be a vector or a structure of vectors')
			end
			notify(obj,'checkData');
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function set.y(obj,value)
			if isstruct(value) && length(value) == 1
				f = fieldnames(value);
				firstLength = [];
				x = [];
				names = {};
				for i = 1:length(f)
					col = value.(f{i});
					if isnumeric(col) && isvector(col)
						if any(size(col)==1) && ~iscolumn(col)
							col = col';
						end
						if isempty(firstLength)
							firstLength = length(col);
						end
						if length(col) == firstLength
							x(:,end+1) = col;
							names{end+1} = f{i};
						end
					end
				end
				if ~isempty(x)
					obj.y = x;
					obj.columnlabels = names;
				end
			elseif isnumeric(value)
				if any(size(value)==1) && ~iscolumn(value)
					value=value';
				end
				value(isnan(value)) = []; %purge nans
				obj.y = value;
			elseif istable(value)
				obj.y = table2array(value);
			end
			
			notify(obj,'checkData');
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function set.columnlabels(obj,value)
			if ischar(value)
				obj.columnlabels = {value};
			elseif iscell(value)
				obj.columnlabels = value;
			end
			
			notify(obj,'checkData');
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function set.cases(obj,value)
			if isempty(value) || (length(value) ~= length(obj.x(:,1)))
				obj.cases = [];
			else
				if size(value,2) > size(value,1)
					value = value';
				end
				for ii = 1:length(value)
					if isempty(value{ii}) || isnan(value{ii})
						value{ii} = 'Unknown';
					end
					value{ii}=regexprep(value{ii},'^\?$','Unknown');
					value{ii}=matlab.lang.makeValidName(value{ii});
				end
				obj.cases = categorical(value);
			end
			notify(obj,'checkData');
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function ret = haveData(obj)
			if isempty(obj.x)
				ret = false;
			else
				ret = true;
			end
        end
        
        % ===================================================================
		%> @brief reset the data to empty
		%>
		%> @param obj this instance object
		% ===================================================================
		function reset(obj)
			obj.x = [];
            obj.y = [];
            obj.xrownames = {};
            obj.yrownames = {};
            obj.cases = [];
            obj.columnlabels = {};
            obj.runStructure = [];
            obj.xdata = [];
            obj.ydata = [];
        end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function doSinglePlots(obj,pn)
			if ~exist('pn','var');pn = obj.pn;end
			wid = 1000;
			hei = 800;
			minp = 0.01;
			maxp = 0.925;
			if obj.isDataEqualLength
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,1).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,2).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,3).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h=figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(2,1).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h=figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(2,2).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h=figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(2,3).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
			else
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,1).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,2).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(1,3).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(2,1).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
				h = figure;
				figpos(1,[wid hei]);
				set(h,'Color',[0.9 0.9 0.9])
				p = copyobj(pn(2,2).axis,h, 'legacy');
				set(p,'Units','Normalized','OuterPosition',[minp minp maxp maxp]);
				
			end
		end
	end
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	methods ( Static = true )
		function [violinplot,medianplot] = violin( x, y, varargin)
			% VIOLIN Violin plot                               
			% Input
			%   x           single numeric value
			%   y           numeric vector
			%
			%   Optional input (key-value pairs, case-insensitive):
			%   'n'         Scalar, number of values of non-parametric kernel
			%               Default = 1000
			%   'facecolor' Color, according to matlab plot color properties
			%               Default = [0 0 0] = black
			%   'facealpha' Transparancy, where 0=fully transparent, 1=fully opaque.
			%               Default = 0.5
			%   'withmdn'   Logical, if true median is plotted as a line. 
			%               Default = false
			%   'linestyle' Linestyle, according to matlab plot linestyle properties
			%               Default = '-'
			%   'linewidth' Linewidth, according to matlab plot linestyle properties
			%               Default = 1
			%   'linecolor' Linecolor, according to matlab plot linestyle properties
			%               Default = 'k'
			%   'style'     1 = opaque outline + tranparent violin patch only
			%               2 = transparent violin patch only
			%               3 = opaque outline only
			%               Default = 1
			%   'side'      'both'  = plot two-sided violin (Default)
			%               'left'  = plot only left-sided violin
			%               'right' = plot only right-sided violin
			%   'rotation'  orientation of violin plots
			%               'vertical'   = violin stands upright
			%               'horizontal' = violin lays down
			%   'scaling'   scalar, scaling of width of violin, in x coordinates.
			%               numerical value
			%               Default = 1
			%   'cutoff'    cutoff density
			%               Default = 1e-3
			%   'kernel'    The type of kernel smoother to use, chosen from 
			%               'normal'  'box', 'triangle', and 'epanechnikov'.
			%               Default = 'normal'
			%   'kernelwidth' scalar, for the bandwidth of the kernel smoothing window.
			%                 The default is optimal for estimating normal densities,
			%               but you may want to choose a smaller value to reveal
			%               features such as multiple modes.
			%
			% EXAMPLE
			% figure(1); clf;
			% rng(17371218);
			% violin(1,randn(1,1000),'KernelWidth', 0.5,'withmdn', 1);
			% violin(0,randn(1,1000) + 6,'KernelWidth', 0.5,'Rotation','horizontal','side','right','facecolor','r','scaling',4);
			%-------------------------------------------------------------------------%              
			for k = 1:2:length(varargin)
				 switch lower(varargin{k})
					  case 'n'
							n = varargin{k+1};
					  case 'facecolor'
							facecolor = varargin{k+1};
					  case 'facealpha'
							facealpha = varargin{k+1};
					  case 'withmdn'
							withmdn = varargin{k+1};
					  case 'linestyle'
							linestyle = varargin{k+1};
					  case 'linewidth'
							linewidth = varargin{k+1};
					  case 'linecolor'
							linecolor = varargin{k+1};
					  case 'style'
							style = varargin{k+1};
					  case 'side'
							side = lower(varargin{k+1});
							if ~ismember(side,{'both','left','right'})
								 error('%s is not a valid value for ''side'' argument in violin',varargin{k});
							end
					  case 'rotation'
							rotation = lower(varargin{k+1});
							if ~ismember(rotation,{'vertical','horizontal'})
								 error('%s is not a valid value for ''rotation'' argument in violin',varargin{k});
							end
					  case 'scaling'
							scaling = varargin{k+1};
					  case 'cutoff'
							cutoff = varargin{k+1};
					  case 'kernel'
							kernel = lower(varargin{k+1});
							if ~ismember(kernel,{'normal','box','triangle','epanechnikov'})
								 error('%s is not a valid value for ''kernel'' argument in violin',varargin{k});
							end
					  case 'kernelwidth'
							kernelwidth = varargin{k+1};
					  otherwise
							error('''%s'' is not a valid key input argument for violin',varargin{k});
				 end
			end

			%% Set default parameter values                     
			if ~exist('n','var') || isempty(n);n = 1000;end
			if ~exist('facecolor','var') || isempty(facecolor)
				 facecolor = 'k';
			elseif isnumeric(facecolor) && length(facecolor)==1
				 facecolor = ones(1,3)*facecolor;
			end
			if ~exist('facealpha','var') || isempty(facealpha);facealpha = 0.5;end
			if ~exist('withmdn','var') || isempty(withmdn);withmdn = false;end
			if ~exist('linestyle','var') || isempty(linestyle); linestyle = '-';end
			if ~exist('linewidth','var') || isempty(linewidth); linewidth = 0.5;end
			if ~exist('linecolor','var') || isempty(linecolor)
				 linecolor = 'k';
			elseif isnumeric(linecolor) && length(linecolor)==1
				 linecolor = ones(1,3)*linecolor;
			end
			if ~exist('style','var') || isempty(style);style = 1;end
			if ~exist('side','var') || isempty(side);side = 'both';end
			if ~exist('rotation','var') || isempty(rotation);rotation = 'vertical';end
			if ~exist('scaling','var') || isempty(scaling);scaling = 1;end
			if ~exist('cutoff','var') || isempty(cutoff);cutoff = 1e-3;end
			if ~exist('kernel','var') || isempty(kernel);kernel = 'normal';end
			if ~exist('kernelwidth','var') || isempty(kernelwidth);kernelwidth = [];end
			switch style
				 case 1
					  % do nothing everything is set with current parameter settings
				 case 2
					  linestyle = 'none';
				 case 3
					  facealpha = 0;
				 otherwise
					  warning('Style can only take numbers 1, 2 or 3');
					  style = 1;
			end
			%% Fit y values (nonparametric kernel-smoothing)    
			% fit distribution
			pd = fitdist(y(:),'kernel','kernel',kernel,'width',kernelwidth);

			% ypdf
			ymin = min(y)-(max(y)-min(y));
			ymax = max(y)+(max(y)-min(y));
			ypdf_tmp = linspace( ymin, ymax, n );

			% xpdf
			xpdf_tmp  = pdf(pd,ypdf_tmp);
			if ischar(scaling) && strcmp(scaling,'off')
				 scaling = max(xpdf_tmp);
			end
			if strcmpi(side,'both'); scaling = scaling / 2; end
			% adjust amplitude to max width
			xpdf_norm = xpdf_tmp ./ ( max(xpdf_tmp) / scaling);

			% length of median line
			if withmdn
				 mdnlength = xpdf_norm(find(ypdf_tmp>nanmedian(y),1,'first')-1);
			end

			%% Trim at cutoff                                   
			include = xpdf_tmp > cutoff;
			xpdf_inc = xpdf_norm(include);
			ypdf_inc = ypdf_tmp(include);

			%% Plotting values                                  
			xpdfR = x + xpdf_inc;
			xpdfL = x - xpdf_inc;

			switch side
				 case 'both'
					  xpdf = [xpdfR(:)' fliplr(xpdfL(:)')];
					  ypdf = [ypdf_inc(:)' fliplr(ypdf_inc(:)')];
				 case 'left'
					  xpdf = fliplr(xpdfL(:)');
					  ypdf = fliplr(ypdf_inc(:)');
					  xpdf = [x xpdf(:)' x];
					  ypdf = [ypdf(1) ypdf(:)' ypdf(end)];
				 case 'right'
					  xpdf = [x xpdfR(:)' x];
					  ypdf = [ypdf_inc(1) ypdf_inc(:)' ypdf_inc(end)];
			end

			switch rotation
				 case 'vertical'
					  % do nothing everything is already set
				 case 'horizontal'
					  xpdfVERT = xpdf;
					  ypdfVERT = ypdf;
					  ypdf = xpdfVERT;
					  xpdf = ypdfVERT;
			end

			%% Plot                                             
			violinplot = patch( xpdf, ypdf, facecolor);

			% beautify
			set( violinplot, 'EdgeColor', linecolor );
			set( violinplot, 'LineStyle', linestyle );
			set( violinplot, 'LineWidth', linewidth );
			set( violinplot, 'FaceAlpha', facealpha );

			% plot median
			if withmdn
				 hold on
				 switch rotation
					  case 'vertical'
							switch side
								 case 'both'
									  medianplot = plot( x+[-1 1].*mdnlength, ones(1,2)*nanmedian(y) );
								 case 'left'
									  medianplot = plot( x+[-1 0].*mdnlength, ones(1,2)*nanmedian(y) );
								 case 'right'
									  medianplot = plot( x+[ 0 1].*mdnlength, ones(1,2)*nanmedian(y) );
							end
					  case 'horizontal'
							switch side
								 case 'both'
									  medianplot = plot( ones(1,2)*nanmedian(y), x+[-1 1].*mdnlength );
								 case 'left'
									  medianplot = plot( ones(1,2)*nanmedian(y), x+[-1 0].*mdnlength );
								 case 'right'
									  medianplot = plot( ones(1,2)*nanmedian(y), x+[ 0 1].*mdnlength );
							end
				 end
				 set(medianplot, 'Color', linecolor );
				 set(medianplot, 'LineWidth', linewidth );
			else
				medianplot = [];
			end
		end
	end
	
	%=======================================================================
	methods ( Access = private ) %-------PRIVATE METHODS-----%
	%=======================================================================
		
		% ===================================================================
		%> @brief Function that runs after the checkData event fires
		%>
		%> @param obj this instance object
		% ===================================================================
		function doCheckData(obj, src, evnt)
			obj.salutation([evnt.EventName ' event'],'Event is running...',true);
            
            if size(obj.x,2) > 100
                error('Too many colums of data, remember rows are data and columns are groups!')
            end
            
			if isempty(obj.y) && ~isempty(obj.x)
				obj.y = zeros(size(obj.x));
            end
            
            %autogenerate group names
			if size(obj.x,2) > length(obj.columnlabels)
				for i = length(obj.columnlabels)+1 : size(obj.x,2)
					obj.columnlabels{i} = ['DataSet' num2str(i)];
				end
            end
            
			if isempty(obj.cases) && ~isempty(obj.uniquecases)
				obj.uniquecases = [];
            end
            
            %autogenerate row names
% 			if isempty(obj.xrownames) || (length(obj.x) ~= length(obj.xrownames))
% 				ntmp = cell(length(obj.x),1);
% 				for i = 1:length(obj.x)
% 					ntmp{i} = ['Obs_' num2str(i)];
% 				end
% 				obj.xrownames = ntmp;
%             end
% 			if isempty(obj.yrownames) || (length(obj.y) ~= length(obj.yrownames))
% 				ntmp = cell(length(obj.y),1);
% 				for i = 1:length(obj.y)
% 					ntmp{i} = ['Obs_' num2str(i)];
% 				end
% 				obj.yrownames = ntmp;
%             end
            
			
% 			if ~isempty(obj.x)
% 				if ~isempty(obj.cases) && obj.isDataEqualLength
% 					obj.xdata = dataset({obj.x,obj.columnlabels{:}},'ObsNames',obj.xrownames);
% 					ntmp = dataset({obj.cases,'Groups'});
% 					obj.xdata = horzcat(obj.xdata, ntmp);
% 				else
% 					obj.xdata = dataset({obj.x,obj.columnlabels{:}},'ObsNames',obj.xrownames);
% 				end
% 			end
% 			if ~isempty(obj.y)
% 				if ~isempty(obj.cases) && obj.isDataEqualLength
% 					obj.ydata = dataset({obj.y,obj.columnlabels{:}},'ObsNames',obj.yrownames);
% 					ntmp = dataset({obj.cases,'Groups'});
% 					obj.ydata = horzcat(obj.ydata, ntmp);
% 				else
% 					obj.ydata = dataset({obj.y,obj.columnlabels{:}},'ObsNames',obj.yrownames);
% 				end
% 			end
			if length(obj.index) > size(obj.x,2)
				obj.index = [];
			end
			
		end
		
		% ===================================================================
		%> @brief
		%>
		%> @param obj this instance object
		% ===================================================================
		function [xcol,ycol,cases,xouttext,youttext] = removeOutliers(obj,xcol,ycol,cases,xouttext,youttext)
			xmean=nanmean(xcol); %initial values before outlier removal
			ymean=nanmean(ycol);
			xstd=nanstd(xcol);
			ystd=nanstd(ycol);
			switch obj.dooutlier
				case 'limit'
					idx1=find(xcol > obj.outlierlimit);
					idx2=find(ycol > obj.outlierlimit);
					if ~isempty(idx1) || ~isempty(idx2)
						if length(xcol)==length(ycol)
							idx=unique([idx1;idx2]);
							xouttext=sprintf('%0.3g ',xcol(idx)');
							youttext=sprintf('%0.3g ',ycol(idx)');
							xcol(idx)=[];
							ycol(idx)=[];
							if ~isempty(cases); cases(idx)=[]; end
						else
							if length(idx1)+2 < length(xcol)
								xouttext=sprintf('%0.3g ',xcol(idx1)');
								xcol(idx1)=[];
							end
							if length(idx2)+2 < length(ycol)
								youttext=sprintf('%0.3g ',ycol(idx2)');
								ycol(idx2)=[];
							end
						end
					end
				case 'quantiles'
					idx1=qoutliers(xcol);
					idx2=qoutliers(ycol);
					if length(xcol)==length(ycol)
						idx=idx1+idx2;
						xouttext=sprintf('%0.3g ',xcol(idx>0)');
						youttext=sprintf('%0.3g ',ycol(idx>0)');
						xcol(idx>0)=[];
						ycol(idx>0)=[];
						cases(idx>0)=[];
					else
						xouttext=sprintf('%0.3g ',xcol(idx1>0)');
						youttext=sprintf('%0.3g ',ycol(idx2>0)');
						xcol(idx1>0)=[];
						ycol(idx2>0)=[];
					end
				case 'rosner'
					idx1=obj.rosnerOutliers(xcol,obj.outlierp,obj.rosnern);
					idx2=obj.rosnerOutliers(ycol,obj.outlierp,obj.rosnern);
					if ~isempty(idx1) || ~isempty(idx2)
						if length(xcol)==length(ycol)
							idx=unique([idx1;idx2]);
							xouttext=sprintf('%0.3g ',xcol(idx)');
							youttext=sprintf('%0.3g ',ycol(idx)');
							xcol(idx)=[];
							ycol(idx)=[];
							if ~isempty(cases); cases(idx)=[]; end
						else
							xouttext=sprintf('%0.3g ',xcol(idx1)');
							youttext=sprintf('%0.3g ',ycol(idx2)');
							xcol(idx1)=[];
							ycol(idx2)=[];
						end
					end
				case 'grubb'
					[~,idx1]=obj.grubbOutliers(xcol,obj.outlierp);
					[~,idx2]=obj.grubbOutliers(ycol,obj.outlierp);
					if length(xcol)==length(ycol)
						idx=unique([idx1;idx2]);
						xouttext=sprintf('%0.3g ',xcol(idx)');
						youttext=sprintf('%0.3g ',ycol(idx)');
						xcol(idx)=[];
						ycol(idx)=[];
						if ~isempty(cases); cases(idx)=[]; end
					else
						xouttext=sprintf('%0.3g ',xcol(idx1)');
						youttext=sprintf('%0.3g ',ycol(idx2)');
						xcol(idx1)=[];
						ycol(idx2)=[];
					end
				case '2SD'
					mtmp=repmat(xmean,length(xcol),1);
					stmp=repmat(xstd,length(xcol),1);
					idx1=abs(xcol-mtmp)>2*stmp;
					
					mtmp=repmat(ymean,length(ycol),1);
					stmp=repmat(ystd,length(ycol),1);
					idx2=abs(ycol-mtmp)>2*stmp;
					if length(xcol)==length(ycol)
						idx=idx1+idx2;
						xouttext=sprintf('%0.3g ',xcol(idx>0)');
						youttext=sprintf('%0.3g ',ycol(idx>0)');
						xcol(idx>0)=[];
						ycol(idx>0)=[];
						if ~isempty(cases); cases(idx>0)=[]; end
					else
						xouttext=sprintf('%0.3g ',xcol(idx1>0)');
						youttext=sprintf('%0.3g ',ycol(idx2>0)');
						xcol(idx1>0)=[];
						ycol(idx2>0)=[];
					end
				case '3SD'
					mtmp=repmat(xmean,length(xcol),1);
					stmp=repmat(xstd,length(xcol),1);
					idx1=abs(xcol-mtmp)>3*stmp;
					
					mtmp=repmat(ymean,length(ycol),1);
					stmp=repmat(ystd,length(ycol),1);
					idx2=abs(ycol-mtmp)>3*stmp;
					if length(xcol)==length(ycol)
						idx=idx1+idx2;
						xouttext=sprintf('%0.3g ',xcol(idx>0));
						youttext=sprintf('%0.3g ',ycol(idx>0));
						xcol(idx>0)=[];
						ycol(idx>0)=[];
						if ~isempty(cases); cases(idx>0)=[]; end
					else
						xouttext=sprintf('%0.3g ',xcol(idx1>0));
						youttext=sprintf('%0.3g ',ycol(idx2>0));
						xcol(idx1>0)=[];
						ycol(idx2>0)=[];
					end
				otherwise
					
			end
		end
		
		% ===================================================================
		%> @brief Converts properties to a structure
		%>
		%>
		%> @param obj this instance object
		%> @return out the structure
		% ===================================================================
		function out=toStructure(obj)
			out = struct();
			fn = fieldnames(obj);
			for j=1:length(fn)
				out.(fn{j}) = obj.(fn{j});
			end
		end
		
		% ===================================================================
		%> @brief Prints messages dependent on verbosity
		%>
		%> Prints messages dependent on verbosity
		%> @param obj this instance object
		%> @param in the calling function
		%> @param message the message that needs printing to command window
		% ===================================================================
		function salutation(obj,in,message,force)
			if ~exist('force','var'); force = false; end
			if obj.verbose==true || force == true
				if ~exist('in','var')
					in = 'undefined';
				end
				if exist('message','var')
					fprintf(['---> getDensity: ' message ' | ' in '\n']);
				else
					fprintf(['---> getDensity: ' in '\n']);
				end
			end
		end
		
		% ===================================================================
		%> @brief Converts properties to a structure
		%>
		%>
		%> @param obj this instance object
		%> @return out the structure
		% ===================================================================
		function [avg,error] = stderr(obj,data,type,onlyerror)
			avg=mean(data, 'omitnan');
			if nargin<4
				onlyerror=0;
			end
			if nargin<3
				type='SE';
			end
			
			if size(type,1)>1
				type=reshape(type,1,size(type,1));
			end
			
			switch(type)
				
				case 'SE'
					err=nanstd(data);
					error=sqrt(err.^2/length(data));
				case '2SE'
					err=nanstd(data);
					error=sqrt((err.^2/length(data)))*2;
				case 'SD'
					error=nanstd(data);
				case '2SD'
					error=(nanstd(data))*2;
				case '3SD'
					error=(nanstd(data))*3;
				case 'V'
					error=nanstd(data).^2;
				case 'F'
					if max(data)==0
						error=0;
					else
						error=nanvar(data)/nanmean(data);
					end
				case 'C'
					if max(data)==0
						error=0;
					else
						error=nanstd(data)/nanmean(data);
					end
				case 'A'
					if max(data)==0
						error=0;
					else
						error=nanvar(diff(data))/(2*nanmean(data));
					end
			end
			
			if onlyerror==1
				avg=error;
			end
		end
		
		% ===================================================================
		%> @brief rosner outlier removale
		%>
		%>
		% ===================================================================
		function index = rosnerOutliers(obj, y, alpha, k)
			index = [];
			y = y(:);
			n = length( y );
			if nargin < 2; alpha = 0.05; end
			if nargin < 3; k = 1; end
			R = zeros( k+1, 1 );
			% sort deviations from the mean
			ybar = mean( y );
			[ ys, is ] = sort( abs( y - ybar ));
			% calculate statistics for up to k outliers
			for i = 0:k
				yy = ys(1:n-i);
				R(i+1) = abs( yy(n-i) - mean(yy) ) / std(yy);
			end;
			% statistical test to find outliers
			index1 = [];
			imax=0;
			for i = 1:k
				%
				pcrit=1-(alpha/((2*(n-i+1))));
				t=tinv(pcrit, n-i-1);
				lambda(i)=(n-i)*t./sqrt(((n-i-1+t^2)*(n-i+1)));
				%
				if R(i) > lambda
					index=is(n-i+1:end);
					index1 = [ index1 is(n-i+1) ];
				end
			end
		end
		
		% ===================================================================
		%> @brief grubb outlier removale
		%>
		%>
		% ===================================================================
		function [b,idx,outliers] = grubbOutliers(a,alpha,rep)
			if nargin == 1
				alpha = 0.05;
				rep = 0;
			elseif nargin == 2
				rep = 0;
			elseif nargin == 3
				if ~ismember(rep,[0 1])
					error('Please enter a 1 or a 0 for optional argument rep.')
				end
			elseif nargin > 3
				error('Requires 1,2, or 3 input arguments.');
			end
			
			if isempty(alpha)
				alpha = 0.05;
			end
			
			b = a;
			b(isinf(a)) = NaN;
			
			%Delete outliers:
			outlier = 1;
			while outlier
				tmp = b(~isnan(b));
				meanval = mean(tmp);
				maxval = tmp(find(abs(tmp-mean(tmp))==max(abs(tmp-mean(tmp)))));
				maxval = maxval(1);
				sdval = std(tmp);
				tn = abs((maxval-meanval)/sdval);
				critval = zcritical(alpha,length(tmp));
				outlier = tn > critval;
				if outlier
					tmp = find(a == maxval);
					b(tmp) = NaN;
				end
			end
			if nargout >= 2
				idx = find(isnan(b));
			end
			if nargout > 2
				outliers = a(idx);
			end
			if ~rep
				b=b(~isnan(b));
			end
			return
			
			function zcrit = zcritical(alpha,n)
				%ZCRIT = ZCRITICAL(ALPHA,N)
				% Computes the critical z value for rejecting outliers (GRUBBS TEST)
				tcrit = tinv(alpha/(2*n),n-2);
				zcrit = (n-1)/sqrt(n)*(sqrt(tcrit^2/(n-2+tcrit^2)));
			end
		end
		
		% ===================================================================
		%> @brief notBoxplot is a boxplot alternative
		%>
		%>
		% ===================================================================
		function varargout=notBoxPlot(obj,y,x,jitter,style)
			% notBoxPlot - Doesn't plot box plots!
			% Rob Campbell - January 2010
			narginchk(1,5)
			if nargin==1
				disp('Wrong number of inputs!')
				return
			end
			
			if isvector(y), y=y(:); end
			
			if nargin<3 || isempty(x)
				x=1:size(y,2);
			end
			
			if nargin<4 || isempty(jitter)
				jitter=0.5; %larger value means greater amplitude jitter
			end
			
			if nargin<5 || isempty(style)
				style='patch'; %Can also be 'line' or 'sdline'
			end
			style=lower(style);
			
			if jitter==0 && strcmp(style,'patch')
				warning('A zero value for jitter means no patch object visible')
			end
			
			if isvector(y) & isvector(x) & length(x)>1
				x=x(:);
				
				if length(x)~=length(y)
					error('length(x) should equal length(y)')
				end
				
				u=unique(x);
				for ii=1:length(u)
					f=find(x==u(ii));
					h(ii)=obj.notBoxPlot(y(f),u(ii),jitter,style);
				end
				
				
				%Make plot look pretty
				if length(u)>1
					xlim([min(u)-1,max(u)+1])
					set(gca,'XTick',u)
				end
				
				if nargout==1
					varargout{1}=h;
				end
				
				return
				
			end
			
			if length(x) ~= size(y,2)
				error('length of x doesn''t match the number of columns in y')
			end
			
			%We're going to render points with the same x value in different
			%colors so we loop through all unique x values and do the plotting
			%with nested functions. No clf in order to give the user more
			%flexibility in combining plot elements.
			hold on
			[uX,a,b]=unique(x);
			h=[];
			for ii=1:length(uX)
				f=find(b==ii);
				h=[h,myPlotter(x(f),y(:,f))];
			end
			hold off
			%Tidy up plot: make it look pretty
			if length(x)>1
				set(gca,'XTick',unique(x))
				xlim([min(x)-1,max(x)+1])
			end
			if nargout==1
				varargout{1}=h;
			end
			function h=myPlotter(X,Y)
				
				SEM=SEM_calc(Y); %Supplied external function
				SD=nanstd(Y);  %Requires the stats toolbox
				mu=nanmean(Y); %Requires the stats toolbox
				
				%The plot colors to use for multiple sets of points on the same x
				%location
				cols=hsv(length(X)+1)*0.5;
				cols(1,:)=0;
				jitScale=jitter*0.9; %To scale the patch by the width of the jitter
				
				for k=1:length(X)
					thisY=Y(:,k);
					thisY=thisY(~isnan(thisY));
					thisX=repmat(X(k),1,length(thisY));
					
					if strcmp(style,'patch')
						h(k).sdPtch=patchMaker(SD(k),[0.85 0.85 0.85]);
					end
					
					if strcmp(style,'patch') || strcmp(style,'sdline')
						h(k).semPtch=patchMaker(SEM(k),[0.8 0.7 0.7]);
						h(k).mu=plot([X(k)-jitScale,X(k)+jitScale],[mu(k),mu(k)],'-r',...
							'linewidth',2);
					end
					
					%Plot jittered raw data
					C=cols(k,:);
					J=(rand(size(thisX))-0.5)*jitter;
					
					
					h(k).data=scatter(thisX+J, thisY, 'o', 'MarkerEdgeColor', C,...
						'MarkerFaceColor', C+(1-C)*0.65,'MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.75);
				end
				
				if strcmp(style,'line') | strcmp(style,'sdline')
					for k=1:length(X)
						%Plot SD
						h(k).sd=plot([X(k),X(k)],[mu(k)-SD(k),mu(k)+SD(k)],...
							'-','color',[0.2,0.2,1],'linewidth',2);
						set(h(k).sd,'ZData',[1,1]*-1)
					end
				end
				
				if strcmp(style,'line')
					for k=1:length(X)
						%Plot mean and SEM
						h(k).mu=plot(X(k),mu(k),'o','color','r',...
							'markerfacecolor','r',...
							'markersize',10);
						
						h(k).sem=plot([X(k),X(k)],[mu(k)-SEM(k),mu(k)+SEM(k)],'-r',...
							'linewidth',2);
						h(k).xAxisLocation=x(k);
					end
				end
				
				function ptch=patchMaker(thisInterval,color)
					l=mu(k)-thisInterval;
					u=mu(k)+thisInterval;
					ptch=patch([X(k)-jitScale, X(k)+jitScale, X(k)+jitScale, X(k)-jitScale],...
						[l,l,u,u], 0);
					set(ptch,'edgecolor','none','facecolor',color)
				end %function patchMaker
				
				function ci = bootCI_calc(vect, nboot, fhandle, alpha)
					mm = fhandle(vect);
					xci = bootci(nboot,{fhandle,vect},'alpha',alpha);
					ci = xci;
					ci = [xci(1) mm xci(2)];
					ci = diff(ci);
					ci = max(ci);
				end
				
				function sem=SEM_calc(vect, CI)
					narginchk(1,2)
					
					if isvector(vect)
						vect=vect(:);
					end
					
					
					if nargin==1
						stdCI = 1.96 ;
					elseif nargin==2
						CI = CI/2 ; %Convert to 2-tail
						stdCI = abs(norminv(CI,0,1)) ;
					end
					
					sem = ( (nanstd(vect)) ./ sqrt(sum(~isnan(vect))) ) * stdCI ;
				end
				
				function tint=tInterval_Calc(vect, CI)
					narginchk(1,2)
					if isvector(vect)
						vect=vect(:);
					end
					if nargin==1
						CI = 0.025; %If no second argument, work out a 2-tailed 5% t-interval
						stdCI=tinv(1-CI, length(vect)-1);
					elseif nargin==2
						CI = CI/2 ; %Convert to 2-tail
						stdCI=tinv(1-CI, length(vect)-1); %Based on the t distribution
					end
					if stdCI==0
						error('Can''t find confidence iterval for 0 standard deviations!')
					end
					tint =  ( (nanstd(vect)) ./ sqrt(sum(~isnan(vect))) ) * stdCI ;
				end
			end %function myPlotter
		end %function notBoxPlot
		
		% ===================================================================
		%> @brief Sets properties from a structure or normal arguments,
		%> ignores invalid properties
		%>
		%> @param args input structure
		%> @param allowedProperties properties possible to set on construction
		% ===================================================================
		function parseArgs(obj, args, allowedProperties)
			
			%lets make allowedProperties from the class properties
			if ~exist('allowedProperties','var') || isempty(allowedProperties)
				ptmp = properties(mfilename);
				otmp = '';
				for i = 1:length(ptmp)
					if isempty(regexpi(ptmp{i},obj.ignoreProperties))
						if i == 1
							otmp = ptmp{i};
						else
							otmp = [otmp '|' ptmp{i}];
						end
					end
				end
				allowedProperties = otmp;
			end
			
			allowedProperties = ['^(' allowedProperties ')$'];
			
			while iscell(args) && length(args) == 1
				args = args{1};
			end
			
			if iscell(args)
				if mod(length(args),2) == 1 % odd
					args = args(1:end-1); %remove last arg
				end
				odd = logical(mod(1:length(args),2));
				even = logical(abs(odd-1));
				args = cell2struct(args(even),args(odd),2);
			end
			
			if isstruct(args)
				fnames = fieldnames(args); %find our argument names
				for i=1:length(fnames)
					if ~isempty(regexp(fnames{i},allowedProperties, 'once')) && isempty(regexp(fnames{i},obj.ignoreProperties, 'once')) %only set if allowed property
						obj.salutation(fnames{i},'Configuring setting in constructor');
						obj.(fnames{i})=args.(fnames{i}); %we set up the properies from the arguments as a structure
					end
				end
			end
			
		end
			
	end
	
end

